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Abstract 

> 

00 I Stochastic chemical systems with diffusion are modeled with a reaction- 

^^ ' diffusion master equation. On a macroscopic level, the governing equation 

^f-s . is a reaction-diffusion equation for the averages of the chemical species. 

,^ I On a mesoscopic level, the master equation for a well stirred chemical 

^^ ' system is combined with Brownian motion in space to obtain the reaction- 

0|0 , diffusion master equation. The space is covered by an unstructured mesh 

and the diffusion coefficients on the mesoscale are obtained from a finite 
element discretization of the Laplace operator on the macroscale. The 
S^ ' resulting method is a flexible hybrid algorithm in that the diffusion can be 

handled either on the meso- or on the macroscale level. The accuracy and 
the efficiency of the method are illustrated in three numerical examples 
inspired by molecular biology. 
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Abbreviations 



Abbreviation 


Expanded form 


mRNA 


messenger ribonucleic acid 


CME 


chemical master equation 


PDF 


probability density function 


SSA 


stochastic simulation algorithm 


RRE 


reaction rate equations 


ODE 


ordinary differential equation 


RDME 


reaction-diffusion master equation 


NRM 


next reaction method 


NSM 


next subvolume method 


BD 


Brownian dynamics 


PDE 


partial differential equation 


RDE 


reaction-diffusion equation 


nD 


n space dimensions 


FEM 


finite element method 


SPDE 


stochastic PDE 


FVM 


finite volume method 



Table 0.1: Abbreviations in the paper in order of appearance. 



1 Introduction 



Intrinsic noise in biochemical networks can have a large impact on the macro- 
scopic behavior of biological cells |lll.l29l.l32l.l36l.l38l|. An example is the regulation 
of the transcription of genes to messenger RNA (mRNA) where genes are present 
in one or two copies and the copy number of mRNA is small. The facts that the 
copy number is a small nonnegative integer and that there is a probability that a 
certain reaction will occur when two molecules meet make a discrete, stochastic 
description of the system necessary. 

The state of the system is the number of molecules of each participating 
species. Usually, the assumption is that the system is well stirred such that there 
is no spatial dependence of the distribution of the chemical species. Then the 
chemical master equation (CME) is the governing equation for the probability 
density function (PDF) of the state of the chemical system 16|, l27|]. A trajectory 
of the biochemical system is simulated by randomly choosing a reaction and 
then u pda ting the state vector in Gillespie's Stochastic Simulation Algorithm 
(SSA) [l9|], further developed in [g, |9|, [l8|, \2f^ to improve the efficiency and to 
handle systems with slow and fast time scales. The concentrations of the species 
at a macroscopic level are often approximated by the reaction rate equations 



(RRE) which is a deterministic system of nonhnear ordinary differential equations 
(ODEs). This approach works well only when the number of molecules is large, 



a condition which is often violated inside living cells [21 . 

There are many biochemical systems where the spatial inhomogeneity of the 
species cannot be neglected. Such systems are no longer well stirred since the 
transport of the molecules through the solvent is slow compared to typical re- 
action times |[28|] or since some reactions are strongly localized. The correlation 
length, i.e. the length scale on which the system can be regarded as s pat i ally 
homogeneous, is now much shorter. Examples are found in [71, M, llO, uA, |33 



where both the stochastic properties and the spatial distribution are necessary 
to explain experimental data. If the diffusion at a molecular level is treated as 
a special set of reactions in the CME, then we arrive at the reaction-diffusion 



master equation (RDME) [16|, Ch. 8], [27|, Ch. XIV], [SOj. This is an equation for 



the time evolution of the PDF of the state of the system in the same manner as 
the CME but the dimensionality of the problem is much higher. 

In the RDME, the geometry is partitioned into computational cells and in 
each cell there is a stochastic variable representing the number of molecules of 
each species. If the number of species is N and the number of cells is K then 
the RDME is an equation for the scalar PDF in KN dimensions and time. For 
example, for a small problem with five species and 100 nodes in a mesh, the 
dimensionality is 500 and a direct solution of the RDME is obviously out of the 
question. The only feasible way is to generate samples from the RDME and 
collect statistics in a Monte Carlo fashion. Examples where the SSA has been 
applied to react ion- diffusion systems in one space dimension can be found in 



[3|, |42|. An efficient version of the SSA is the next reaction method (NRM) [18|. 
An implementation of the NRM, specially developed for diffusive systems, is the 
next suhvolume method (NSM) [lO(] , implemented in the freely available computer 



software MesoRD 23|. 



More accurate modeling may be necessary if the number of molecules is very 
low. In Brownian dynamics (BD), the separate paths of single molecules are 
tracked and they may react with other molecules in their vicinity [l|, |39|, |47|. The 



reaction and diffusion of the particles are simulated using the Green's function 



of the Smoluchowski and diffusion equation in [47j. The BD approach is possible 
only if the total number of molecules is small and becomes inefficient when there 
are many nonreactive collisions for each reactive one. The RDME and BD are 
compared in [T]]. 

The corresponding macroscopic equation for the concentrations of the species 
is a partial differential equation (PDE): the react ion- diffusion equation (RDE). 
This is the equation solved in applications such as combustion [37] and the model 
is appropriate when the number of molecules is large and stochastic effects can 
be neglected. 

In this paper, we develop a method for simulating the RDME on an unstruc- 
tured mesh consisting of triangles in two space dimensions (2D) or tetrahedra in 



3D. Unstructured meshes have the advantage of approximating curved inner and 
outer boundaries much more accurately than Cartesian meshes do. The diffusion 
coefficients at the meso level are chosen to be consistent with the discretization 
with the finite element method (FEM) on the mesh converging to the diffusion 
operator at the macro level when the cell size vanishes. With a proper com- 



putational mesh 17j, the FEM approximation yields positive probabilities for a 



particle to jump into the adjacent cell. The time integration of the RDME is split 



according to Strang |41| into two parts in a hybrid method: the diffusion and the 
chemical reactions. First, a macroscopic diffusive step is taken for a subset of 
the species with the diffusion operator at the macro level using the unstructured 
primal mesh. Then the stochastic reactions and the stochastic diffusion for the 
remaining species are advanced in the dual cells of the mesh with SSA at the 
meso level. If all species in all cells are treated at the meso level, then we exactly 
sample the RDME. 

Hybrid methods for a well stirred s yst em with a mesoscopic-macroscopic ap- 



proximation are found in [23, 12J]. In |34|. an efficient Monte Carlo method for 
the diffusion equation is described. In [26|, the mesoscopic diffusion coefficients 
are derived from the discretization of the Laplacian on a 2D Cartesian mesh with 
an interior boundary. With unstructured meshes, we propose in this paper to 
obtain those coefficients from a proper FEM discretization. 

The outline of the paper is as follows. In Section [21 the RDME is stated 
and the relation to the RDE is discussed. The diffusion coefficients are derived 
from the discrete Laplacian in Section [31 The first and second moments of the 
distribution of molecules in a system without chemical reactions are derived in 
Section [H The hybrid method coupling the meso- and macroscales and the op- 
erator splitting in time are found in Section [3 Three examples in 2D are found 
in the section with numerical results. It is shown that the suggested mesoscopic 
diffusion yields trajectories converging to the macroscopic diffusion equation. It 
is also shown that the hybrid method we propose accurately samples the RDME 
at a fraction of the time needed for a full simulation. Conclusions are drawn in 
the final section. 



2 Reaction-Diffusion Master Equation 

Assume that the computational domain Q in space is partitioned into computa- 
tional cells Cj, j = 1, . . . ,K, such that the cells do not overlap and they cover 
the whole domain 

d nCj =^,iy^ j, and uf^^ Cj = Q. 

Furthermore, assume that there are A^ chemically active species Xij, i = 1, . . . ,N, 
in the K cells, j = 1, . . . ,K. The state of the system is the array x with N x K 
components Xij. The jth column of x is denoted by x.j and the ith row by Xj.. 



The non-negative integer Xij is thus the copy number of species i in cell j. The 
time dependent state is changed by chemical reactions occurring between the 
molecules in the same cell and by diffusion where molecules move to adjacent 
cells. In the reactions, the species interact vertically in the array x and in the 
diffusion, the interaction is horizontal. The RDME governs the time evolution of 
the PDF p, where p(x, t) is the probability to be in state x at time t. 

2.1 Chemical reactions 

A reaction r in a cell j is a transition from one state x.j before the reaction to the 
state x.j = x.j — Ur after the reaction. The state-change vector n^ of a reaction 
is a vector with small integer numbers of length A^ independent of j. There is 
a reaction probability or propensity Wr that reaction r will take place in a cell 
depending on the state x.j. A chemical reaction in cell j can be written 

x.j y x.j, x.j = x.j + 11^. l^-lj 

In a system without diffusion, the PDF for the molecular distribution in Cj 
satisfies the CME. Let n.^ be split into two parts 

rir = n+ + n^, n^^ = max{nri, 0), n'^ = mm{nri, 0). 

Then the CME for p is, see 0, Ch. 7], [27, Ch. V], 

9p(x, t) 



dt 



= Mp{^,t)= (2.2) 

K R 

^ ^ Wr{x.j + nr)p{x.i, . . . , X.j + n,,, . . . , X.x, t) 

j=l r=l 

x.j+n,, >0 

K R 



J^ U7^(x.j)p(x,t), 



j=l r=l 

x.j— n^ >0 

where the constraints on x are defined elementwise. These constraints are intro- 
duced in order to avoid unfeasible reactions and will be dropped in the following 
as is customary. 

A simple reaction in Ck is 

Xik + Xjk — > Xik, wi(x.fc) = cikXikXjk-, (2.3) 

where conventionally, we use uppercase letters to denote molecule names, while 
lowercases are used for counting the number of molecules of a certain species. The 
transition vector rii is zero except for three components: Un = n^ = l,nii = 
—1. The propensity Wi has the same form for all cells k, whereas the reaction 
coefficient Ci^ scales with the area or volume \Ck\ of the cell such that cik = Ci/\Ck\, 
where ci is a constant. 



2.2 Diffusion 

Suppose now that there are no chemical reactions but only diffusion in the system. 
Then the mesoscale model of the diffusion of species i from one cell Ck to another 
cell Cj can be written as a chemical reaction (cf. (12. 3p ) 

Xik -^ Xij, Vkj{yii.) = qkjXik- (2.4) 

It is understood that qkj is non-zero only for those cells that are connected and 
Qjj = 0. The form of the propensity Vkj and the diffusion coefficient q^j are the 
same for all species here but q^j may depend on i and can be different for small 
and large molecules. We can write 

(Ikj = 7t^, (2-5) 

where 7 is the macroscopic diffusion constant, hkj is a measure of the local length- 
scale and qkj is dimensionless but still depends on the precise shapes of the cells 
Ck and Cj. The interpretation of qkj as the inverse of the expected value of the 
first exit time for a single molecule from cell Ck to Cj makes it clear why no 
simple formula exists except for very regular cells. The molecular movement is 
often modeled by the Ito-diffusion 

d^ = adWt, (2.6) 

where ^(t) is the position and Wt is a Wiener-process. The relation between 7 
in (12.51) and a in (12.61) is then simply 7 = a'^/2. 

With this notation, the master equation for the diffusion in fl2.4p can be 
written in the same manner as the CME in (Q, see [l6|, Ch. 8], [27, Ch. XIV], 

„ . . N K K 

op[:s., t) sr^ \-^ \-^ / \ / \ 

— ^7 — = Z^Z^ / ykji^i- + mfcj)p(xi., . . . , Xj. -I- nifcj, . . . , xat., t) 

1=1 k=i j=i 

-ffej(xi.)p(x,t). (2.7) 

The corresponding transition vector rrifcj is zero except for two components: 
mkj,k = 1 and rrikjj = -1. 

By combining (12.21) and (12.71) . we arrive at the RDME for a chemical system 
with reactions and diffusion 

^P^ = Mpi^,t)+Vpi^,t). (2.8) 

at 

We now give a few comments on the validity of the RDME. Denote the molec- 
ular reaction radius by pji and the shortest average life time of the molecular 
species by Tmm |2|]. Then the requirement for the size /i of a cell is 

Pr < /i^ < a7T-min, (2.9) 



where a. is O {X) and depends on the cell shape and the dimension |lOl . l27l |. 
Firstly, the upper bound guarantees that the mixing in a cell by diffusion is 
sufficiently fast for the molecules to be homogeneously distributed there. Another 
interpretation is that with slow diffusion, a better spatial resolution is necessary 
since the solution becomes less smooth. Secondly, there is also a lower bound 
on h for the modeling at a meso level to be meaningful. This lower bound 
guarantees that association and disassociation events can be properly localized 
within a computational cell. The model breaks down if we let /i — > despite 
the fact that there is often a meaningful stochastic partial differential equation 
(SPDE) in the limit. This SPDE, however, only remains valid as the first term 



in a system size expansion [16|, Ch. 8.2]. It is understood with such an expansion 
that we are working on a scale where the cell size appears small but is still large 
enough to contain many molecules. 

Let Cfc be a cell with a part of its boundary on d^. Since there is no transport 
of molecules out from f2 at C^ with the reactions in (12.41) . the corresponding 
boundary condition at the macro level is a Neumann condition. Most results for 
the numerical solution of diffusive systems are derived under Dirichlet boundary 
conditions and we therefore now discuss the associated conditions for the RDME. 
Suppose that fi is surrounded by a reservoir such that the number of molecules is 
fixed at the boundary d^. Then x.^ in those cells is given by the reservoir data. 
If Xij = 0, then Vjk = in (12. 4p and there is no diffusion of molecules from the 
boundary to the interior cells, (cf. (12. 4p ) 

Xik -^ Xij = 0, Vkj{xi.) = QkjXik- (2.10) 

The transition vector nikj is zero except for rukj^k = 1- 

2.3 Relation to macroscopic equations 

Define the concentrations (ptj of species i in cell Cj at a macroscopic level as the 
expected values of \Cj\~^Xij. The RRE for (pij is obtained from the CME (12.21) by 



deriving equations for the mean values of Xij as in |16l. l27l]. The system of ODEs 
defining the RRE is 

R 



^^ u;,{<f>.^)^-J2T^^Wr{\C,\ct^.,). (2.11) 

If there are no reactions but only diffusion, then a similar set of macroscopic 
equations may be derived. From the similarity between (12.21) and (12.71) we have 

K 

fe=i """' 






Yl w^\^kj(pik - I y^qjk I (pij- (2.12) 

fc=i \^^\ 




The diffusion equation (12.121) has the form 

'-§^;D^l (2.13) 

for each i (cf. (12.51) ). The diffusion matrix has the elements 7Djfc = Q'fej|Cfc|/|Cj|, j ^ 
k, and -fDj^ = - J2k^j qjk- 

The positive coefficients q^j in our diffusion model are non-zero only if cell k 
and j share a common point in ID, a common edge in 2D, or a common facet 
in 3D. The diffusion matrix D therefore has a sparsity pattern matching the 
connectivity of the partitioning of Q into computational cells. 

Assume for now that the diffusion is isotropic on a Cartesian lattice in ID 
with constant cell size h so that the probability qkj to move from cell Ck to cell Cj 
is equal to the probability qjk to move in the opposite direction. With isotropy, 
D has a particularly simple structure and is symmetric. For example, if fi = [0, 1] 
with h = 1/K and diffusion q we have that 

4>ii = q{-4>ii + 4>i2), 

4>ij = q{(j)i,j-i - 2(j)i, + 0,,,+i), J = 2, . . . , K - 1, (2.14) 

4>iK =q{4>i,K-i-4>iK)- 

Let q = '-f/h? as in (12. 5p and let h —>■ 0. Then the solution of (12.141) converges to 
the solution (f)i{x,t) of the diffusion equation in ID 

If boundary values 0ji and (j)iK are given, then the solution converges to the 
solution of the PDE in (12.151) with boundary conditions 0j(O,t) = (pn = Qio and 
0i(l,t) = (t)iK = gn for some Qio^gn > 0. 

Suppose that the interior of il = [0, 1] x [0, 1] is covered by square cells of size 
h X h in 2D or that Q = [0, 1] x [0, 1] x [0, 1] in 3D and partitioned into cubic 
cells of size h x h x h. In both cases, D is symmetric and will approximate the 
Laplacian A with Neumann boundary conditions. With the normal derivative of 
(pi at the boundary dQ written as dcpi/dn, the solution 0j converges when /i ^ 
to the solution of 

^ = 7A0, in fi, ^ = on dU, (2.16) 

ot on 

If (pi is given data in the boundary cells, then the boundary conditions in (12.161) 
will be Dirichlet type, (pi = Qi on dQ. 

The macroscopic approximation of the diffusive part of (12. 8p with a vanishing 
cell size in a Cartesian mesh thus satisfies a diffusion equation (I2.15p . The macro- 
scopic counterpart to Ai in (12.80 is the RRE (12.110 . A macroscopic concentration 
(pi affected by both chemical reactions and diffusion fulfills a RDE 

^ = 07,(0) + 7A0„ t = l,...,N. (2.17) 



A Cartesian mesh for discretization of A0 has many advantages but is im- 
practical for curved inner and outer boundaries of Q. A RDME on unstructured 
meshes is proposed in Section [3l 

2.4 Relation to microscopic equations 

If we consider pure diffusion in ID, it is possible to directly compare the coeffi- 
cients obtained from the finite element discretization with the expected value of 
the first exit time of Brownian motion from a finite interval. Using linear basis 
functions in a FEM discretization of (12.151) on a mesh with vertices Xj and cell 
sizes hj = Xj — Xj_i, in the interior of Q the non-zero entries of the stiffness 
matrix 5* and the mass matrix M are 



I ,. I ,. 1 



i J''j+i ^ "-i "-i+i (2.18) 



After mass lumping the coefficients corresponding to jumps from node j to its 
neighbors in (12.51) are given by 

27 27 

Consequently, the exponentially distributed waiting time for the next event at 
node j has the expected value 

27 
?ij-i + ?ij+i = 7-r • (2-20) 

On a uniform grid we recover the jump coefficients 7//^^ and the parameter 2'y/h'^ 



used in 23 



Interpreted as the average time for a Brownian particle starting at node Xj to 
reach either of its neighbors, we can compare the value in (I2.20p with the actual 
expected value of the first exit time r = inf{t : ^t ^ (xj_i,Xj+i)} where C,t is 
defined in (12. 6p . A straightforward application of Dynkin's formula (48|, Ch. 7.4] 
yields 

^^.H-^ = ^, (2.21) 

in accordance to (12.201) . The probabilities to exit at Xj^i and Xj+i respectively are 
given by hj+i/{hj + /ij+i) and hj/{hj + /ij+i), and using this we recover the jump 
coefficients f l2.19p . In this sense, the coefficients of the mesoscale model obtained 
from the dicretization of the macroscale equation is consistent with the microscale 
description. It is worth noticing however, that r is not generally exponentially 
distributed 0, p. 212]. 



3 Diffusion coefficients 




Figure 3.1: The primal mesh (thick hnes) with the vertices (small circles), the 
dual mesh (thin lines) and the bisectors of the triangles (dashed lines). 

Consider a part of an unstructured mesh in 2D covering Q with a polygonal 
boundary dQ in Figure [XT! The primal mesh consists of triangles with the vertices 
in the corners. The cells Ck in the dual mesh are polygons and in the interior 
of Q, the center of Ck is the vertex k. The edges of the polygon coincide with a 
part of the bisectors of the triangles or with the boundary dQ. The corners of 
an inner Ck are the barycenters of the triangles and the midpoints of the edges 
from its center vertex. In ID, the cell in the primal mesh is a line segment with 
a vertex in both ends and the dual mesh also consists of line segments shifted 
with respect to the primal mesh and a vertex in the center. The primal cell is a 
tetrahedron and the dual cell is a polyhedron in 3D. 

Similarly, the dual cells Ck in a Cartesian structured mesh in 2D with vertices 
{ih,jh),i = 0, . . . ,K, j = 0, . . . ,K, are the cells in the staggered mesh. In the 
interior of Q, the vertices are in the center of Ck- 

A finite element discretization of the RDE in (12.171) with continuous, piece- 
wise linear basis and test functions on the primal unstructured mesh generates 
a system of equations to solve for the nodal values (f). The components of 4> are 



4>ij{t), the concentration of species i in vertex j at time t 4J, Ch. 15] but can 



also be interpreted as the mean value of the concentration of species i in the dual 
cell Cj. The system of equations is [44, Ch. 14] 



M(/) = f (0) + 750. 



(3.1) 



The mass matrix M and the stiffness matrix 5* are symmetric, M is positive 
definite and S is negative semi-definite with Neumann boundary conditions and 
negative definite with Dirichlet conditions. The reaction terms are represented 



10 



by the nonlinear term f . If /i is a suitable measure of the sizes of the triangles 
in the mesh, then the error in the FEM solution of the Dirichlet problem in the 
L2-norm is O {h^) 0, Ch. 14]. 

Order the unknowns in according to a linear index so that 

0=(0i.,---,0,.,---,0^.f- (3.2) 

Then M and S are block diagonal matrices where the blocks M and S are identi- 
cal, small mass and stiffness matrices. The system (13.11) is simplified by introduc- 
ing mass lumping of M and f . Let A be a diagonal matrix with diagonal blocks 
A with Ajj = Ylk=i ^ik and let similarly f be the lumped version of f . In ID, 
Ajj is the length of the dual cell with vertex j in the center, see (12.181) . Similarly 
in 2D, Ajj is the area \Cj\ of the dual cell in Figur eJXTl with vertex j in the center 



4J, Ch. 15] and in 3D, Ajj is the volume of Cj 451]. The simplified system (13.11) 



IS now 



(i) = Cj{ci))+-ibct), u = A'^^,b = A-^S, (3.3) 

which is our approximation of (12.17^ . The solution </> of (13. 3p is also generally 
second order accurate in 2D [351 . 

Let D denote a block on the diagonal of D. Then D = A^^S and an off- 
diagonal component Dj^. is non-zero only if two vertices j and k are connected 
by an edge. With Neumann conditions, the diagonal blocks D satisfy Djj = 
— '^k^jDjk so that J2k=i^jk = fo^ every j. In other words, the constant 
vector ei = (1, 1, ... , 1)^ is in the null-space of D and is the right eigenvector 
with eigenvalue zero. The corresponding left eigenvector 62 has the diagonal 
elements of A as components, e2j = Ajj. Let 



Xij 



^Xijp{x,t) = \Cj\ 



denote the expected value of the number of molecules of species i in cell j. In a 
system without reactions, by (13.31) 

d d ^ d ^ d ^ 

j=i j=i j=i 

i.e. the total number of molecules of each species are conserved by the diffusion. 
This is not the case with Dirichlet conditions, where Djj < — X]a;=^7 ^jk and D is 
non- singular. 

By the calculation in Section 12.41 in ID, we find that Dj^ >Oiik=j — 1 
or k = j + 1 and zero for the other non-diagonal entries. If the mesh in 2D is 



a Delaunay triangulation [17[, then the Voronoi cells are close to the cells in the 



dual mesh defined by Figure [3TT1 Two triangles share an edge in Figure [3^ . The 
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(a) 



(b) 



Figure 3.2: The critical angles in 2D (a) and 3D (b) for positive off-diagonal 
elements in the stiffness matrix. 

sum of the two opposing angles a and /3 in the triangles is less than or equal to vr 
in a Delaunay triangulation and Dj^ > when j ^ k 46|]. In 3D, assume that the 
dihedral angle at an edge between two facets of a tetrahedron, see Figure 13.2b . 
is non-obtuse {a < 7r/2) for all tetrahedra in the mesh. Then the elements 



of D all have the right sign [46||. With these properties of the discretization 
and without chemical reactions, a discrete maximum principle is fulfilled 46| for 
Dirichlet boundary conditions ensuring that with non-negative initial conditions 
the solution remains non-negative. The diffusion matrix D generated by the FEM 
discretization of A in fl2.17p has the same properties as the diffusion matrix in 
fl2.13p . A review of these triangulations in three and higher dimensions is found 
in [5|]. They may not be trivial to generate in higher dimensions than 2. 

Discretization matrices of second order PDEs are frequently M -matrices |2C 



10.3]. Henceforth, we shall assume the following slightly weaker property: 

Assumption 3.1 The diffusion matrix D for Dirichlet and Neumann boundary 
conditions fulfills for j ^ k, 



D 



jk 



> 0, Djj < 0, 



K 

fc=i 



D,k < 0. 



The last inequality is an equality for Neumann conditions. 



The macroscopic elements in 'jD in (13. 3p define the coefficients q^j in the 
mesoscopic model of diffusion (12.41) . The diffusion matrix D in fl2.13p has the 
form 'jD = A~^QA, where Qjk = qkj when j ^ k. Since D in (12.131) and (13.31) 
are identical we have 

Q = ^SA-^ = -fD^. (3.5) 

Except for the case when all cells have the same size, Q is generally unsymmetric. 



12 



The concentrations (pij defined as the expected values oi Xij/Ajj with the PDF 
in (12. 8p will satisfy (13. 3p . The molecules at the meso level can jump between dual 
cells with a point (ID), edge (2D), or facet (3D) in common since Djk > there. 
The connectivity graph of D tells in which direction a molecule can diffuse and 
the positive elements of D are inversely proportional to the expected value of the 
first exit time to leave the cell, cf. (12.51) and Section [2^ 

An alternative to the finite element discretization of the RDE in (I2.17P is to 
use the finite volume method (FVM). Here, the averages of the concentrations 
(f)ij in the dual cells are the degrees of freedom. Then 



^ij 



dt 



IMic, Ci mi^c, (3g) 

= -- / cui((/)) dV+-— 7n ■ V0i dS, 

mi Jcj mi JdCj 

where Gauss' theorem has been used and fi is the normal of dCj. The reac- 
tion term in Cj is approximated by uji{4>.j). The gradient V0i is needed on the 
boundary of Cj and different approximations are possible. A simple one is to 
let fi ■ V0i ~ {(pik — 4>ij)/hjk where hj^ is the distance between vertex j and a 
neighboring vertex k. The resulting diffusion term is in Cj, 

/ ^ Cjkifpik — <Pij)/hjk, 

k 

with summation over those vertices with an edge between j and k and where Cjk 
equals the size of the edges or facets of dCj adjacent to the same edge. Then 
the coefficients in this discretization are interpreted as the diffusion coefficients 
in (12. 7p in the same manner as in the FEM case. The difference between (13.60 
and (13. 3p lies in the approximation of the diffusion term. 

Convergence of the FVM to the analytical solution is proved for certain dis- 



cretizations of the gradient in [13|, |25| but the quality of the approximation seems 
to depend critically on the quality of the mesh [43|. This is one of the rea- 
sons why we prefer the FEM approach. Another reason is that FVM is perhaps 
more suitable for problems dominated by convection while chemical systems from 
molecular biology tend to be of diffusive character. 

4 Moments of the diffusion 

The purpose of this section is to prove that the diffusion can be accurately evolved 
deterministically when the number of diffusing molecules is sufficiently large. 
Consider the equations for the moments of x in a system with diffusion and 
Neumann boundary conditions but without chemical reactions. The expected 



13 



values and the covariance matrices satisfy systems of ODEs. These equations are 



derived in [12|, |27|| for general propensities, but with 



being linear (cf. (12 ■4p ). they have a particularly simple structure. 

Applying the formulas in 12|, [l5| or invoking (12.130 . the first moment of the 



number of molecules of species z in a cell is given by 

x^ = Q^l- (4.1) 

This equation is exact since the diffusion propensities v^^ are linear in x and no 
coupling to higher order moments exists. 

The second moments or covariances of any species i between cells j and k are 
denoted by Cj^. The equation for Cj^ is 

K K 

Cjk = Yl <5fc;C'j7 + Yl <5j'Ch + F,k, j,k = l,...,K, (4.2) 

1=1 1=1 

with the driving term F defined by 

K K 

^ik = y^ y^^a/3j"^a/3,fcfa/3(x^.). 
a=l (3=1 

This can be written in matrix form since C is symmetric, 

C = CQ^ + QC + F. (4.3) 

Using the properties of the diffusion propensities, the elements of F are 

K 

^n = XI Qji^ii + Qij^ij^ Fjk = -{QjkXik + QkjXij), j ^ k, (4.4) 
i=i,¥j 

The covariance equation is also exact for the same reason as before. 
The solution of (14.11) can be written 

Xi.(t)=x°.exp(Q^t), (4.5) 

where x° is the initial value at t = 0. The eigenvalues of Q and D are all negative 
except for one which is zero. The corresponding left eigenvector of Q is ei and 
the right eigenvector is 62, both of them defined in Section [3l Hence, 

5ci.{t) = Kie^ + 6yii.{t). (4.6) 
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with an upper bound on ||(5xj.(t)|| given by csexp{\2t) where A2 is the negative 
eigenvalue of Q with smallest magnitude and the norm is the ^2-norm. Therefore, 

lim Xj.(t) = Kje^. (4.7) 

By (13.41) . Xj.ei is constant and we obtain 

K K 

Kieler = x°.ei ^ k^ = J] 4/ ^ Ajj. (4.8) 

The explicit solution of (14.31) is 

C(t) = exp(Qt)Co exp(g'^t) + / exp(Q(t - s))Fexp(g'^(t - s)) ds, (4.9) 

Jo 

where Cq is the initial value of the covariance at t = 0. Using (14.41) . (14.61) . (14. 7p . 
and (13.51) we find that when t — >■ cxo 

Fjj = 7Ki E£i,/^i AiSji/Aii + AjjSij/Ajj = 27Ki Ylf=i,i^j Sji = -2-^KiSjj, 
Fjk = -^Ki{AkkSjk/Akk + AjjSkj/Ajj) = -2-fKiSjk, j ^ k, 

and hence that, 

F = -2-fKiS + 6F, (4.10) 

where \\SF\\ is bounded by c^i? exp(A2t). A bound on the covariance matrix for 
finite time is derived in the next proposition. 

Proposition 4.1 Suppose that Cq = 0. Then 

\\Cm<CF^^^^^^ t\\^4ds, (4.11) 

^^^j^jj Jo 

for some bounding constant cp such that \\F\\ < Ci?||xj. ||. 

Proof. Let S = A-^/^SA-^/\ Then 

exp(gt) = exp{^SA-H) = A^/^ exp{^St)A-^/^. 

The symmetric matrix 76* has the same eigenvalues as Q. Let the unitary ma- 
trix U have the eigenvectors of S as columns and let A have the corresponding 
eigenvalues on the diagonal. A bound on the exponential of Qt with r > is 



exp 



{Qr)\\ < \\A^^^\\\\A-^/^\\Uexp{-fAT)U^ 



< II exp(7Ar) || max \/ A^^j min \/ Aj^ < max yAjj/ min \/Ajj . 
j j j j 
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The same bound is valid for exp{Q'^T). Hence, 

/ exp(Q(t-s))Fexp(Q^(t-s)) rfs 
Jo 

I . . , ^ max, /ii,- 



< / F niaxAi-i/mmA,'! (is < cp , / Xj. Lis. 



D 



Remark. Using fl4.9p and (I4.10p . one can show that ||C(t)|| is bounded when 
t-^oo. If Co = 0, then \\C{t)\\ ~ Ki for large t. U 

It follows from the proposition that the variance of Xj. is proportional to ||xj.|| 
in a bounded time interval so that the standard deviation is proportional to 
a/||xj.||. When t is large then Xj. in (14.71) is of the same order as Kj, a weighted 
average of the initial x°. (14.81) . The standard deviation is proportional to ^/E'i 
according to the remark after Proposition 14.11 Therefore, the quotient between 
the standard deviation and the expected value is small for large copy numbers 
implying that the expected value is indeed a good approximation. Conversely, 
if the number of molecules in a cell is small, then a description in terms of 
expectation values should not be used. 

5 Time integration and hybrid diffusion 

For a discretization parameterized by the cell size /i, the time to compute a 
trajectory of a system with SSA spent in the diffusion part of (12.81) is proportional 
to x7//i^, where x is the total number of diffusing molecules (cf. (12.41) and (12. 5p ). 
It follows that diffusion is the only event with a total intensity that increases 
with increasing spatial resolution. In order to avoid this, we propose to split the 
diffusion operator T) into two parts with the species with low copy numbers in 
Vl and the species with high copy numbers in Vh- The diffusion in Vh can then 
be advanced in time macroscopically. 

Order the species Xj such that Xj, « = !,..., X^, have low copy numbers and 
Xj, z = Ni + 1, . . . , X, have high numbers and let 

K K 

6 = ^ ^Wfcj(xj. + mfcj)p(xi., . . . , Xj. + nifcj, . . . , Xat., t) - Vkj{^i.)p{^, t), 

k=l 3=1 

Nl 

VLp{i^,t) = Y^i,,VH = V-VL. (5.1) 

i=l 

The operator on the right hand side of (12.81) can be written 

^P^ = [M+ Pi]p(x, t) + Vhp{^, t). (5.2) 
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It follows from Section H] that the effect of diffusion on the species with many 
molecules in each cell is well approximated by the mean-field equations. In the 
numerical solution procedure, the second part of ( 15.20 is therefore first advanced 
half a time step At/2 with the macroscopic diffusion. Then the first part is 
integrated a full step At and finally the macroscopic diffusion is applied for half 



a time step again. This is the Strang splitting procedure [4l| to solve (15.21) and 
below we give conservative conditions under which both steps preserve the non- 
negativity of the solution. 

For many relevant cases, a single trajectory gives sufficient insight into the 
stochastic reaction-diffusion system but the same procedure also works well for si- 
multaneous simulation of M trajectories. An application could be to approximate 
the PDF by following an ensemble of trajectories with state vectors x™(t), m = 
1, . . . , M. Then p is reconstructed according to 

M 



m=l ^ 

In order to advance the trajectories At in time from t" to t""*"^, each trajectory 
is first integrated in time from t" to t" -|- At/2 by solving (14. ip for the species 
i = Nl + 1, . . . , A^. The time derivative in (14.11) is discretized by the trapezoidal 
(or Crank-Nicolson) method of second order temporal accuracy and the new state 
x" for each trajectory is the solution of 

(/ - \MQ){^t"Y ={!+ ^Atg)(x,.(t"))^, z = AT^ + 1, . . . , AT. (5.4) 
Alternatively, a scheme with an error O (At) is the Euler backward method 

{I-\AtQ){^T''r={Mnf- (5.5) 

The time step in the Strang splitting must be sufficiently small to resolve the 
shortest time scale of the reactions Tmin in (12. 9p . Thus, we can take At ~ Tmin 
and by ([23]) and [ij 



h^ < a7At. (5.6) 

For such a time step, an explicit method for integration of (14. ip would be very 
inefficient. Only an implicit method suitable for stiff problems, such as the trape- 
zoidal method, can efficiently advance the solution of (14. ip in time in a stable 
manner. 



The next step is to evolve the M trajectories with SSA [19| for one time step 
At using the reduced master equation 

^^^ = A<p(x, t) + Vlp{^, t). (5.7) 
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If At is short then there may be no events in many of the reahzations of the 
process and some computational work will be wasted. 

The final step is to evolve each trajectory x™' half a time step again using 
the macroscopic diffusion. The result is an approximative sample of p in (15. 2p at 
time r+\ 

The error due to the Strang splitting and the time discretization by the trape- 
zoidal method is O (At^). Similarly, the error due to the FEM approximation of 
the diffusion is O (h^) 0, Ch. 7]. 

It is more difficult to estimate the error induced from using the macroscopic 
diffusion without being overly pessimistic. A bound on the local single trajectory 
error can be obtained as follows. Write ||C(At)|| = (9 (At||xj.||) by Proposi- 
tion 14. 1[ Then the local stochastic error (standard deviation) relative to Xj. is 
O (At^/^||xj.||~^/^) for i > Nl since only species with high copy numbers partic- 
ipate in the macroscopic diffusion. This simple bound, however, gives no global 
estimate since reasonable and sufficient stability properties of the system are dif- 
ficult to prescribe. When using deterministic diffusion for some species it must 
simply be regarded as given a priori that the diffusion noise for those variables 
has little or no impact on the system as a whole. Note that computing averages in 
order to approximate expected values will generally make this error substantially 
smaller since the macroscopic diffusion is exact in expectation (cf. (14.11) ). 

In a system with only diffusion, the algorithm has the following two properties. 

Proposition 5.1 Assume that D satisfies Assumption \3.1[ that x'^(O) > for 
all trajectories m = 1, . . . , M, and that At < ^min/67 for the trapezoidal method 
in ^5.4\ ), where /imin is the minimal distance between a vertex and the opposing 
edge in a triangle in the mesh. Then in a system without chemical reactions, the 
copy numbers in the trajectories computed by the hybrid algorithm remain non- 
negative for t > 0. For the Euler backward method 1(5. 5]) . there is no time step 
restriction for non-negativity. 

Proof. If D satisfies Assumption 13. H then Q satisfies the assumption by (13. 5p . 
Let (xj.(t"))"^ = A^/^y" in a symmetrization of (15.41) . Then 

(/ - 5)y"+^/^ = (/ + S)y" = g, ^ = Q.h-iAtA-^I^SA-^'^. 

The symmetric matrix S has components Sjk > 0, j 7^ /c, and Sjj < and 



/ — S" is positive definite since S is negative semi-definite. By [4jy, Lemma 15.4], 



(/ - S)^^ > and by ^ Theorem 15.6] if At < h'ijG-f and y" > 0, then the 
right hand side g is non-negative. Consequently, 

y"+V4^(j_5)-lg>0, 

and therefore x" > 0. The only differences for the Euler backward method 
are that the right hand side g equals y" and is always non-negative and S is 



twice as large. The intermediate SSA-step also preserves the non-negativity as 
do the final step. Thus, the copy numbers of the species computed by the hybrid 
algorithm remain non-negative. D 

Remark. The upper bound on At in the proposition is quite restrictive considering 
the requirements for resolving the reactive time scale in the splitting in (15. 6p . In 
practice the solutions in the macroscopic diffusion step stay non-negative with 
much longer At since the mean values are large for the species involved in VhP- 

D 

Proposition 5.2 In a system with only diffusion and 'Yj^. Djk = 0, j = 1, . . . , i^, 
the total number of molecules of each species in a trajectory is constant. 

Note that the exact solution to the equations for the concentrations has the 
same property in (13. 4p . 

Proof. The vector ei satisfies -Dei = and e^Q = ejD^ = 0. In the first step of 
the hybrid algorithm (15.41) . we have 

ef (/ - lAtQ)i^r'^r = e^i^^'^r = e^I + ^Atg)(x,.(t"))^ = ef (x,.(t"))^. 
Consequently, 



M K M K 

^n+1/4 ^~~^ ^~~^ ^m, 71+1/4 



m=l j=l m=l j=l 



and Sj, the sum of the copy number of species i over all cells, is thus preserved 
by (15. 4p . This sum over all cells is also preserved by diffusion simulated by SSA 
in every trajectory in the intermediate step in the hybrid algorithm. Finally, Si 
is preserved in the last step of the Strang splitting. D 

The accuracy and efficiency of the algorithm are evaluated in the next section. 



6 Numerical results 

The algorithm for the RDME in Section O is applied to the diffusion equation 
and to two different systems from molecular biology. The convergence of samples 
from the mesoscopic diffusion model to the solution of the macroscopic equation is 
illustrated in Section [Ol The method is applied to a model of a bi-stable reaction 
network in Section 16.21 Finally, in Section 16.31 we illustrate the potential of the 
hybrid method by comparing it to a purely stochastic simulation. The meshes. 



the stiffness and the mass matrices are generated using the PD&toolbox [3l[ in 
MATLAB. 
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Figure 6.1: The triangular mesh with 123 vertices (left) and the logarithm of 
the error vs. the number of vertices (right). The ^2-norm of the error (solid) is 
compared to the asymptotic rate of convergence (dashed). 



6.1 Diffusion 

The diffusion equation with Neumann boundary conditions and initial data 

ut = l{u^x + Uyy), m.n= [-0.5, 0.5] x [-0.5, 0.5], 7 = 10~^ 

du 

— = 0, on dVL, u{x, y, 0) = 100(1 - cos(27rx)), 

has the analytical solution 

Uaix, y, t) = 100(1 - cos(27rx) exp{-A-f7r^t)). (6.1) 

The solution Ud is computed with a FEM discretization of the space deriva- 
tives with mass lumping and integrated in time by the trapezoidal method as 
in (15. 4p . The error Ud — Ua in the vertices in the ^2-iiorm is O (At^) + O {h"^), 
see Section [5l The behavior of the spatial error at t = 1 is confirmed in Fig- 
ure 16. 1[ Asymptotically, when the number of vertices increases the theoretical 
rate is obtained. 

A stochastic simulation of the diffusion with m molecules in the mesh is 
compared with the analytical and the FEM solutions in Table 16. 1[ Two meshes 
are used: one with 33 vertices and a maximum length h^ax of an edge of a triangle 
equal to 0.5 and one with 123 vertices and hmax = 0.25 (see Figure EH]). Each 
stochastic simulation starts with 100 molecules distributed according to u{x, y, 0). 
The average concentrations Um in the cells Ck are computed for M trajectories 
such that m = lOOM. Then Um is compared to Ua and Ud in the vertices at t = 1. 
The weighted vector norms in £2 and i^o are defined by 



\u 



Zl^'i'^ii' 



\u\ 



max \u 



'i\ 



(6.2) 



The differences in these norms divided by the system size 100, Sa = (um — Ua) /lOO 
and Sd = (um—Ud) /lOO, for two different discretizations are collected in Table lOl 
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^ 


2 




^CXD 




i^max 


= 0.5 


iT'max - 


= 0.25 


'^rnax 


= 0.5 


'T'max 


= 0.25 


logio ^ 


5a 


5d 


Sa 


5d 


Sa 


5d 


Sa 


5d 


2 


.038 


.038 


.1 


.1 


.52 


.53 


4.5 


4.5 


3 


.019 


.018 


.029 


.029 


.32 


.31 


1.2 


1.2 


4 


.0053 


.0053 


.0073 


.0073 


.091 


.096 


.23 


.23 


5 


.0018 


.0016 


.0029 


.0030 


.031 


.017 


.16 


.17 


6 


.0011 


.00057 


.0011 


.0010 


.016 


.0093 


.047 


.042 


7 


.00083 


.00015 


.00039 


.00027 


.013 


.0023 


.014 


.010 



Table 6.1: The relative difference between the stochastic solution and the analytic 
solution 6a or the FEM solution 5^ for different mesh sizes hmax and total number 
of molecules m. 

The error is expected to behave as O (/i^aa;) + ^ (^~^^^) and this is what we 
observe in the table. The difference between Um and Ud decays with increasing 
m at the correct rate in both norms. In the example in Figure 16.11 the £2-error 
is 8.7 ■ 10~^ when hmax = 0.5 and 2.8 • 10~^ when hmax = 0.25 explaining the 
difference between 6a and 6d in Table 16.11 When m is large then the dominant 
term in 6a is the discretization error. 



6.2 Domain separation in a bi-stable system 

In this section we simulate a model of a bi-stable system, previously investigated 



using the freely available software MesoRD [23|] in [10|. The model consists of 
eight chemical species participating in twelve reactions, see Table [6^ Being based 
on a double negative feedback mechanism, in the spatially homogeneous case the 
system switches between states with mostly A molecules and states where B is 
dominating. The model is used to illustrate and explain the observation that 
global bi-stability can be lost in a spatially dependent system due to domain 
separation, when the diffusion is slow. 



Ea^Ea + A Eb^Eb + B E_ 



'^A 



B ^ EaB Eb + A^ EbA 



EaB -\- B ^ EaB2 EbA ^ EbA2 

kd kd 



A 



k4 



B 



k^ 



Table 6.2: The chemical reactions of the bi-stable model. The constants take the 
values fci = l50s-\ ka = 1.2 x lO^s-^M-^, k^ = Ws'^ and k^ = 6s-\ 



We have used an implementation of the NSM [lOj, with support for unstruc- 
tured meshes added by us. The code is written in C, and wrapped in a MATLAB 
mex-file. This approach makes the definition of the geometry, the meshing and 
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the matrix assembly convenient. 









(a) t = Os 



(b) t = 2.08s 



(c) t = 8.32s 



Figure 6.2: Snapshots of the time evolution of one trajectory with diffusion coef- 
ficient 7 = 2 X lO^^^m'^/s. The concentrations of the species A and B are found 
in the upper and lower rows, respectively. Dark areas indicate regions with a 
higher concentration. Patches with the system in different phases are formed at 
t = 8.32s. 



The unstructured mesh has K = 8849 nodes, giving a minimal dual cell area 
of 8.13 X lO^^^m^. The mesh quality is not perfect; a few off diagonal elements 
in the stiffness matrix fail to be non- negative, giving slightly wrong diffusion 
rates locally. This, however, does not seem to have any profound impact on the 
behavior of the simulated system when compared to simulations on structured 
triangular meshes where all coefficients are of the correct sign. The boundaries are 
reflecting, corresponding to a Neumann boundary condition in the finite element 
formulation of the macroscopic equation. 

The time evolution of the system is simulated on a circle with radius 3 x 10~^m 
in Figure ES] with 7 = 2 x lO'^^m^/s. There is a domain separation with many 
A molecules in the lower right part and many B molecules in the upper left part. 

In Figure 16.31 the same system is simulated with fast diffusion, 7 = 1 x 
10"^^r7i^/s, and in this case the system does not separate into domains with 
different phases. At the end time of this simulation, the system is in a state where 
A molecules dominate, and the system behaves much like in the homogeneous 
case. 
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(a) t = Os 



(b) t = 2.16s 



(c) t = 7.62s 



Figure 6.3: The same system as in Figure [6^ is simulated with 7 = 1x10 ^^m^/s. 
Here the system has not separated into domains with different phases. 

6.3 The hybrid method — metaboUtes and enzymes 

A model of a biochemical network with two metabolites A and B and two enzymes 
Ea and Eb from [40] is simulated in this example. The domain fi is a disc with 
radius pmax = vr~^/^ ^ 0.564 and area |f2| = 1. The reactions are summarized 
in Table 16.31 Initially, the concentrations a and b are constant and the enzyme 
concentrations, ca and e^, are zero in every cell. 



Wi[a,eA} 
W2{h,eB) 



A 



5^0 A+B 



k2 



E, 






'^B 



W-i 



W4 



kBCB/il + hkj^) w^ib) 



■ H {0.2 ~p)kE J (l + ak],'] 
H{p-OA)kEjil + bk^') 



Table 6.3: Reaction channels for the network. The concentrations of the species 
A, B, Ea, and Eb in a cell are a, b, ca, ^b- The reaction constants are kA = kB = 
3C, p = 0.002C, k2 = 0.0005C^ fc^^ = ks^ = 0.5C, kj = 60/C, kn = 30/C and the 
diffusion constant is 7 = 10~^. The domain is covered by i^ = 80 dual cells and 
the scaling with the average size of a cell ( = \^\/K is done to define the unit 
scale of the problem. H denotes the Heaviside function; the enzyme Ea is thus 
produced only in the center of fi, < p < 0.2, and Eb is created only close to 
the boundary, 0.4 < p < Pmax- 
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The chemical reactions and the diffusion of the enzymes are simulated with 
SSA and the diffusion of A and B is modeled by the diffusion PDE in a straightfor- 
ward MATLAB implementation of the hybrid method. The effect of the enzymes 
in different parts of Q is demonstrated in Figure [63] with M = 10^ realizations. 






Figure 6.4: The concentrations of the metabolites A (upper left) and B (lower 
left) and enzymes Ea (upper right) and Eb (lower right) at t = 200. 



A stochastic reference solution 0*, i = 1,2,3,4, averaged over M = 10^ 
trajectories is integrated directly to t = 200 with SSA. Hybrid solutions 0j. *(200) 
are computed for different At with M = 10^ for comparison with (pl. Table [631 
displays the differences between the stochastic and the hybrid solutions. 



At 


0.1 


1 


5 


20 


40 


100 


St 


0.024 


0.024 


0.024 


0.024 


0.025 


0.030 



Table 6.4: The difference between the stochastic solution and the hybrid solution 
in the £2-norm, 6t = maxj 1101. — 0^. *||2/(maXj (/)^j — imuj (p'lj) at t = 200 for 
different At. The stochastic errors dominate in 6t for At < 40. 

The computational work for the stochastic and the deterministic parts of the 
algorithm is compared in Figure 16.51 (a) at t = 200 with 7 = 10~^ and M = 10^. 
The work in the deterministic part is less than the work in the stochastic part 
when At > 1, and decreases as At~^ when At increases. The work for the 
stochastic part tends to a limit since the extra effort for restarting SSA in each 
step becomes negligible when the time step is longer. 

Figure [^751 (b) displays the total CPU-time for stochastic and hybrid solutions 
for different diffusion coefficients 7 at t = 10 with M = 10^. We use At = 5 for 
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the hybrid algorithm, i.e. a time step well below the upper limit for which the 
stochastic errors dominate for 7 = 10~^. The stochastic solutions are integrated 
directly to t = 10. When the diffusion of the molecules is the major part of 
the computational work for larger 7 > 10~^, then replacing the diffusion for 
the species with large copy numbers by the macroscopic diffusion reduces the 
CPU-time by up to 1000 while retaining small differences between the solutions, 
see Table 16.41 The number of molecules of the enzymes rif. = (X3. + X4.)ei is 
about 0.001 of the number of metabolite molecules n^ = (xi. + X2.)ei at t = 10. 
This difference explains the remarkable improvement in efficiency in diffusion 
dominated regimes. From the figure and the discussion, the CPU-times Tssa 
and Thyb for SSA and the hybrid method, respectively, are approximately 

Tssa ~ c^o + Cs7(^m + n^), Thyb ~ cto + CsjUe, 

where Cso and Cho are small and are mainly due to the chemical reactions. The 
concentration of the enzymes increases when the integration is continued to t = 
200. The speedup by the hybrid method with At = 5 is then about 35. Since 
nm/ne ~ 30 this is in agreement with the estimates of Tssa and T^yb above. 



e e — e e 




10" 10' 

time step 




10"' 10"' 

diffusion coefficient 



(b) 



Figure 6.5: Run times versus the time step for the stochastic and deterministic 
part of the hybrid algorithm (a) and run times versus the diffusion coefficient 7 
for SSA and the hybrid algorithm (b). 



7 Conclusions 

The RDME is discretized on an unstructured mesh for better geometric flexibility 
than a Cartesian mesh. The diffusion coefficients in the RDME are derived from 
a FEM discretization of the Laplacian. Stochastic simulation can then be used to 
determine a number of trajectories of the mesoscopic system. If the copy number 
is large for some chemical species, then a hybrid method integrating the diffusion 
part deterministically reduces the computing time substantially, especially when 
the diffusion constant is large. 
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The method is apphed to three different systems. The convergence of a sys- 
tem without reactions to the macroscopic solution is shown in the first example. 
In the second example, we consider a biochemical system with both diffusion and 
reactions. We illustrate that our method can be efficiently applied to a previ- 
ously studied bi-stable system, and that the results obtained with our code are 
in line with those computed on structured meshes with the freely available soft- 



ware MesoRD [23|. In a final example, the performance of the hybrid method is 
compared to SSA for a system with four species. The hybrid algorithm is up to 
three orders of magnitude faster. 

In this paper we have considered examples in two space dimensions only. Re- 
alistic modeling of the reaction networks in e.g. bacteria typically require 3D 
simulations. With our approach, the extension to 3D is straightforward and will 
be reported in a forthcoming paper. The method relies on the FEM discretization 
of the diffusion equation and a variety of existing software can be used to spec- 
ify the geometry, construct the mesh and postprocess the result. Presently, we 
handle diffusion with a uniform diffusion constant but there is no complication in 
considering space-dependent diffusion or adding convection or to let the diffusion 
be different for different species. 

In the spatially homogeneous case, stiffness arises from the presence of fast 
chemical reactions. For the RDME, the stiffness in addition increases with the 
resolution of the mesh. We have proposed a hybrid method to reduce simulation 
time when some species are present in large copy numbers. For the homogeneous 
case, many approximative schemes have been developed for systems with time 
scale separation. For spatially dependent systems, no sharp separation in slow 
and fast events can generally be made. Instead we rather have a continuum of 
scales and multiscale simulation techniques have to be developed anew. 
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